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Results of numerical simulations of a recently derived most 
general dissipative-dispersive PDE describing evolution of a 
film flowing down an inclined plane are presented. They in- 
dicate that a novel complex type of spatiotemporal patterns 
can exist for strange attractors of nonequilibrium systems. It 
is suggested that real-life experiments satisfying the validity 
conditions of the theory are possible: the required sufficiently 
viscous liquids are readily available. 

PACS numbers: 05.45. +b, 47.20.Ky, 03.40Gc 



The phenomenon of pattern formation in nonequilib- 
rium dissipative systems is currently a topic of active 
experimental and theoretical research (see e.g. 0] for a 
recent progress review). Here we report our theoretical 
studies and numerical simulations of a two-dimensional 
(2D) evolution PDE approximating a flow down an in- 
clined plane; it exhibits self-organization of a remarkably 
complex spatiotemporal pattern which then persists in- 
defmetly in this dissipative-dispersive system. In cer- 
tain cases discussed below, such a pattern consists of two 
subpatterns of two-dimensionally (2D) localized surface 
structures. One of these subpatterns is an essentially ID 
arrangement of larger-amplitude bulges on the film sur- 
face which are nearly equidistantly aligned on (a num- 
ber of) straight-line segments; those are surrounded by 
smaller-amplitude bumps, which constitute the second, 
lattice-like subpattern filling up essentially the entire flow 
domain. Each of the two subpatterns moves as a whole; 
their velocities are different, and every elementary struc- 
ture (a bulge as well as a bump) periodically collides with 
those of the other kind. In the collision of a bump with 
a bulge (or with a pair of neighboring bulges), the two 
structures pass through each other similar to the well- 
known ID Korteweg-deVries solitons, returning to their 
pre-collisional shapes and speeds after the interaction. 

Studies of wavy film flows on solid surfaces (the 
"Kapitza problem") have a considerable history. How- 
ever, the nonlinear dynamics of wavy films is far from be- 
ing fully understood (see e.g. j|; see ||] and [Q for recent 
progress reviews). Fortunately, the nonlinear couplcd- 
PDE Navier- Stokes (NS) problem, additionally compli- 
cated with a free boundary, can be reduced to simpler ap- 
proximate descriptions of the wave dynamics for certain 
domains of the parameter space. In the most favorable 
cases, such a description reduces to a single partial dif- 
ferential equation (PDE) governing the evolution of film 



thickness. Recently, we (see ||) have derived the most 
general evolution equation (EE) capable of all-time-valid 
description of a wavy liquid film (of a constant density 
p, kinematic viscosity v, surface tension a, and average 
thickness ho) flowing down an inclined plane. Its dimen- 
sionless form is 



r/t + 4rjr] z + -Sr) z 



2 2 

-cot 9ri yy + -WV\ 



2V 2 r, z = 0. 
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Here r\ is the deviation of film thickness from its average 
value of 1, z and y are the streamwise and spanwise co- 
ordinates, and V 2 = d 2 /dz 2 + d 2 /dy 2 . We have defined 
S = (4i?/5 — cot 0), where is the angle of inclination of 
the plane and R = h§XJ jv is the Reynolds number. Here 
U = gh\ sin 0/(2v) where g is the gravity acceleration, 
and in Eq. (1) W — a/(2pvU) is the Weber number. All 
the dimensionless variables are measured in units based 
on ho, U , and p. Equation (1) describes the film evolu- 
tion in a reference frame moving with the velocity 2U in 
the streamwise direction. 

In this paper, we limit ourselves to the case of 9 = 7r/2, 
i.e. flow down a vertical wall. Then, we can transform 
the EE to a "canonical" form which will contain only one 
control parameter — by rescaling rj = Nrj, z = Lz, y = 
Ly, and t = Tt, where N = 2R/(5W), L = y/5W/(4R), 
and T = (5 3/2 /16)(W/ R) 3/2 . Dropping the tildes in the 
notations of variables, the resulting canonical form of the 
EE is 



+ VVz + VV + e (i] zz + V 4 ?7) = 0. 



The control parameter in this equation is 



= (l/3)VW/5. 



(2) 



(3) 



Equations (1) and (2) have been derived directly from 
the fundamental NS equations by using an iterative pro- 
cedure which is a variation of the so-called multiparam- 
eter perturbation approach (see e.g. Iq], ||, and refer- 
ences therein; an earlier, more limited application of mul- 
tiple independent perturbation parameters appears e.g. 
in fl). In addition to its leading to the most general 
EE, another advantage of this technique is that it yields 
the least restrictive conditions of theory validity. For the 
present case, they require that the following two dimen- 
sionless parameters be independently small: 
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a = y/RfW < 1 and j3 = ^R^/W < 1. (4) 

From the linear stability theory, the third (third-order 
derivative) term of Eq. (^|) is purely dispersive, while all 
other linear terms are dissipative. Different limiting cases 
of Eq. (|2|) reproduce some known nonlinear equations, 
such as the 2D version of Korteweg-de Vries (KdV) equa- 
tion for e — > and 2D version of Kuramoto-Sivashinsky 
(KS) equation for e — * oo (see also The ID limit 

(d y = 0) of Eq. (2) is essentially the well-studied (see 
e.g. references in [gj) equation first introduced by Kawa- 
hara M. 

To exhibit interesting spatial behavior, a system should 
be sufficiently "large" . For the periodicity domain < 
y < 2np, < z < 2nq in our simulations of Eq. (||) , we 
chose 5 < p < 16 and 16 < q < 80. We used spatial grids 
of up to 32p x 32q nodes, with the Fourier pseudospec- 
tral method for spatial derivatives and with appropri- 
ate dealiasing. Time marching was done (in the Fourier 
space) by using Adams- Bashforth and/or Runge-Kutta 
methods. We checked the results by refining the space 
grids and time steps; by verifying the volume conserva- 
tion, J ndydz — 0; etc. A typical simulation run took 
~ 10 6 time-steps. 

The initial values of rj were chosen independently at 
each node from the interval [—0.05, 0.05] with uniform 
probability distribution. Due to the dissipativity of Eq. 
(|J), the system evolves to an attractor, and so essentially 
"forgets" the initial conditions. For large values of the 
"dissipativity" parameter e ^> 1, as soon as the flow ap- 
proaches its asymptotic state, the observed film surface 
is irregular in space and time; no spatial patterns are 
evident. The chaotic character of the attractor is indi- 
cated by the positive largest Liapunov exponent [which 
we found, similar to Deissler [flO||, by numerically inte- 
grating, along with Eq. (2), the linear equation that 
governs small disturbances of the solution of Eq. (2)]. 
This is in accordance with the fact that in the limit 
e — > oo, Eq. (2) reduces to a 2D generalization of the 
Kuramoto-Sivashinsky equation, whose solutions on ex- 
tended spatial domains are known to exhibit chaotic at- 
tractors. Regarding the transient behavior on the way 
to the attractor, our simulations of Eq. (2) with small e 
corresponding to the parameter values of the experiment 
have shown agreement |3j with their transient "3D" 
patterns and pattern transitions, including checkerboard 
patterns, synchronous instabilities, and solitary waves. 

The main focus of the present communication is the 
presence of highly nontrivial orderly patterns in time- 
asymptotic states for the strongly dispersive cases, e -C 1. 
Figure 1 shows snapshots of the film surface at large times 
for three different sets of parameter values. (The fact 
that by those times the systems have approached their 
asymptotic states is clear, e.g., from the corresponding 
plots of the evolution of "energy" J rj 2 dydz [see Fig. 2 
corresponding to Fig. 1(a)]. We will speak of such nu- 



merically identified time-asymptotic states as attractors, 
although one needs to be cautious here: it is known that 
such extended systems may sometimes exhibit long tran- 
sients. We find the largest Liapunov exponent to be pos- 
itive in this case as well, suggesting a strange attractor.) 

There are two subpatterns in Fig. 1(a): The V-shaped 
formation consisting of 13 large-amplitude bulges aligned 
into two straight lines moves as a whole downstream with 
a certain velocity, and the small-amplitude lattice-like 
subpattern of bumps moves uniformly as well, but in 
the opposite (in our reference frame) direction. [Simi- 
lar segregation of coherent structures into two subpat- 
terns of different amplitudes is also seen for the non- 
square, large-aspect-ratio domains, Figs. 1(b) and 1(c).] 
This collision-course movement is evident in the cross- 
sectional space-time portrait shown in Fig. 3. Even 
though each bump changes its shape in irregular man- 
ner, the bump maintains its identity. In particular, the 
bumps do not seem to coalesce or break up, and just 
weakly interact with one another. Also, the height of a 
bulge irregularly fluctuates, by an amount which is ap- 
proximately equal to the amplitude of the bumps. 

As a bump runs into a bulge, the bulge's amplitude in- 
creases momentarily, and then decreases again as a bump 
separates from the opposite side of the bulge (see Fig. 3). 
These interactions, unlike the irreversible coalescences of 
1-D pulses — discovered in jll| for a highly nonlinear dis- 
sipative equation — appear to be (almost) reversible, like 
the well-known interactions of ID KdV solitons. 

We note that bulge formations similar to that of Fig. 
1(a) were discovered in |7j for = 25 and p = 64/(27r) 
(they postulated an equation of the form (2) based on an 
equation of the form (1) derived in Ref. E3] for a partial 
case of an inclined film; in fact, that derivation was not 
valid for the vertical film) . But the authors of seem to 
have overlooked (perhaps, because of inadequacy of the 
graphics tools they used) the second, bump subpattern — 
and thus the entire complex, dynamical character of the 
two-component order. 

It is natural to inquire as to how the various quantities 
of the pattern scale with e. We varied e _1 between 25 
and 305 for p = q = 16. In one set of simulations, e~ l 
was gradually decreased from 50 in relatively small steps 
of 5 (to allow the system to "adiabatically" adjust to the 
new parameter value), up to e _1 = 25 — at which point 
the line formations of bulges break down. In another set 
of simulations, e _1 was increased from 50 in steps of 10 
or 15 up to e _1 = 305. In all cases, we find that the 
characteristic width of the bulge as well as the bump is 
of the order of (~) 1 independent of e. The amplitude 
of bulges is also constant, ~ 1, as are the velocity of 
bulges and that of bumps. Only the amplitude of bumps 
changes; it scales like ~ e. 

The V-shaped formation of bulges retains its form 
when e is changed from 1/30 to 1/305. However, the (ab- 
solute value of the) angle <p of each bulge line with the 
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streamwise axis decreases with e, probably approaching 
some asymptotic value in the limit e — > (see Fig. 4; 
since there is no parameters remaining in this limit, the 
asymptotic angle should be just 0). It might be possible 
to explain this dependence by a theory of pairwise inter- 
action of bulges through their (non-axisymmetric) "tails" 
(similar to the theory G3] for 2D chemical- wave spirals). 

When e <C 1, the dissipative terms in Eq. (||) can be 
treated as perturbations ~ e of the 2D KdV equation 

r)t + Wh + V 2 r/ 2 = 0. (5) 

This equation does not seem to have any analytical so- 
lutions. However, by transforming to a reference frame 
moving with a velocity c > [replace rj t with (— cr) z ) in 
the equation], Petviashvili and Yan'kov (m) numerically 
obtained a stationary axially-symmetric solitary- wave so- 
lution. By balancing the first term with the nonlinear 
term, cr\ z ~ Tft) z , and the latter term with the dispersive 
term, the characteristic amplitude and velocity of these 
solutions are found to be r] ~ c and c ~ 1/L 2 where L s 
is the characteristic lengthscale, which is not uniquely 
determined by the KdV equation (^|). However, the two 
dissipative terms of Eq. (||) will change L s on a slow 
time scale, until they balance each other (this essential 
role of the small dissipative terms was revealed in Ref. 
H for the ID case). This selects the soliton of L s ~ 1, 
which results in c ~ 1 and n ~ 1 as well, independent of e. 
These estimates are clearly consistent with the numerical 
results for bulges reported above. 

Motivated by the discovery of the second, small- 
amplitude subpattern, we examined the possibility of 
a corresponding second travelling-wave solution. If we 
transfer to the frame moving with a negative velocity 
c = —a 2 , where a is a (real) constant, there are such 
solutions — with the nonlinear term being as small as the 
dissipative ones. Indeed, the leading-order equation then 
is V 2 r] z + a 2 rj z = 0, which is the well-known Helmholtz 
equation for r\ z . There are solutions oc sin Jy sin K z 
(J 2 +K 2 = a 2 ). The balance between the (small) dissipa- 
tive terms again determines K <~ c ~ 1, and the balance 
of the dissipative terms with the nonlinear term yields 
77 ~ e. We see that these lengthscale, amplitude and 
velocity (including its sign) agree with those observed 
for the bumps in the numerical experiments as described 
above. Note that our assumption of the negative veloc- 
ity is essential: with a positive velocity, one arrives at 
the modified Helmholtz equation, which does not have 
any oscillating solutions. There are only exponential so- 
lutions, which are unsuitable here. [We note that the 
Helmholtz equation has axially-symmetric solutions as 
well, oc Jo(ar) where J is the Bessel function (r is the 
radial coordinate) . This solution is only weakly localized: 
it decays at spatial infinity like a power rather than ex- 
ponentially. There is no such localized solutions in the 
ID case, d y = 0.] 



One would naturally like to find some known types of 
patterns to which those reported here can be compared. 
There are several known cases (see e.g. Ref. [|l5| and 
references therein) of indefinitely long coexistence of dif- 
ferent types of patterns. However, in those cases each 
of the coexisting patterns is confined to its own spatial 
region: its constituent structures do not penetrate inside 
any "alien" pattern region. In contrast, we have seen that 
the bumps constantly pass through the region of bulges. 
Another possibility would be to look at the bulges and 
bumps as two traveling waves. However, in contrast to 
usual cases, the bulge "wave" is confined to an essentially 
ID region, and there is a constant nonlinear interaction 
with the wave of bumps. 

Similar to Eq. (||), we have derived an equation for 
a film flowing down a vertical cylinder (see Ref. |j) and 
references therein). In particular, one can see that if 
the (dimensionless) radius b of the cylinder is not too 
small (6 ^> /3 _1 ), the flow is well approximated by the 
planar-film equation (|l|) [with periodic BC in the az- 
imuthal direction; we note that this also justifies our use 
of spanwise-periodic BCs in the numerical simulations. 
As to the streamwise BC, we believe the solution be- 
comes essentially insensitive to their specific type in the 
limit of large aspect ratio q/p, as e.g. in Fig. 1. It would 
be interesting to check this with spatial-evolution simula- 
tions, such as those already conducted [p| for a different 
EE that coincides with the non-dispersive limit of the 
above-mentioned EE 0]. One finds that with ho ~ 1 
mm, the cylinder (dimensional) radius b ~ 1 cm, and un- 
der parametric conditions a < 1, ,8 < 1, and e <C 1, for 
the waves (evolving as they propagate from the entrance 
end of not-too-long a cylinder to its exit end) to have 
enough time for approaching the attractor stage, the liq- 
uid should be several hundred times as viscous as water. 
For example, it could be glycerin with an admixture of 
water. 

As a general conclusion, numerical 2D simulations of a 
realistic evolution PDE signal that nonequilibrium dis- 
sipative systems can spontaneously form nonperiodic, 
but nevertheless highly ordered spatial patterns (of com- 
pactly localized, soliton-like structures) which are of a 
remarkable complexity. In particular, the novel pat- 
terns consist of subpatterns — each of a different ampli- 
tude and each moving as a whole with its own veloc- 
ity, penetrating through one another. Thus, these sub- 
patterns coexist as the constituent components of the 
overall — microscopically nonsteady but macroscopically 
permanent — self-organized spatiotemporal order charac- 
teristic of the system motion on such an attractor. 

The particular dissipative-dispersive PDEs (1) and 
(2) have been consistently derived from the full Navier- 
Stokes problem to provide a controllably close approx- 
imation to the evolution of a liquid film flowing down 
a vertical plane. The unconventional perturbative ap- 
proach used in this derivation has the advantage of yield- 
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ing the least restrictive conditions of the validity of the 
theory. To satisfy those validity conditions for a possible 
(terrestrial) experiment designed to observe patterns of 
the novel type on a film flowing down a vertical cylinder, 
the film liquid should be much more viscous than water; 
fortunately, suitable liquids are readily available. 
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FIG. 2. Evolution of the surface deviation "energy" 
J rfdydz from an initial small-amplitude "white-noise" sur- 
face to an attractor of Eq. (2). The snapshot Fig. 1(a) was 
taken near the end of this run. Note that the time unit here 
is 50 times that of Eq. (2). 

FIG. 3. Time-sequence of instantaneous surface profiles 
in a fixed vertical cross-section normal to the film (for 
p = q = 16, e _1 = 50; the time shown as is in fact 1.6 x 10 
counting from the start of the run). In particular, it is evident 
that the (large-amplitude) bulges move in one direction and 
(small-amplitude) bumps in the opposite direction. 



FIG. 4. Angle ip between (each) line of bulges and the 
streamwise direction varies with e (p = q — 16) . 

* Internet address: afrenkel@gp.as.ua.edu 
[1] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys 65, 
851 (1993). 

[2] J. Liu, J. B. Schneider, and J. P. Gollub, Phys. Fluids 7, 
55 (1995). 

[3] A. L. Frenkel and K. Indireshkumar, in Mathemati- 
cal modeling and simulation in hydrodynamic stability, 
edited by D.N. Riahi (World Scientific, Singapore, 1996), 
pp. 35-81. 

[4] H. -C. Chang and E. A. Demekhin, Advances in Applied 

Mechanics 32, 1 (1996). 
[5] A. L. Frenkel, Europhys. Lett. 18, 583 (1992). 
[6] D. J. Benney, J. Math. Phys. 45, 150 (1966). 
[7] S. Toh, H. Iwasaki, and T. Kawahara, Phys. Rev. A 40, 

5472 (1989). 

[8] N. J. Balmforth, Ann. Rev. Fluid Mech. 27, 335 (1995). 
[9] T. Kawahara, Phys. Rev. Lett. 51, 381 (1983). 
[10] R. J. Deissler, J. Stat. Phys. 54, 1459 (1989). 
[11] V. I. Kerchman and A. L. Frenkel, Theoret. Comput. 

Fluid Dynamics 6, 235 (1994). 
[12] J. Topper and T. Kawahara, J. Phys. Soc. Japan 44, 663 
(1978). 

[13] I. S. Aranson, L. Kramer, and A. Weber, Phys. Rev. E 

47, 3231 (1993). 
[14] V. I. Petviashvili and V. V. Yan'kov, Dokl. Akad. Nauk. 

SSSR 267, 825 (1982). 
[15] K. Lerman, G. Ahlers, and D. S. Cannell, Phys. Rev. E 

53, R2041 (1996). 
[16] R. J. Deissler, A. Oron, and Y. C. Lee, Phys. Rev. A 43, 

4558 (1991). 



FIG. 1. Snapshots of the time-asymptotic film surface 
self-organized in simulations of Eq. ^, for three different cases 
(bulges move down the page here; for convenience of presen- 
tation, different axes may have different scales; in reality, all 
"bulges" and "bumps" have small slopes and are nearly ax- 
isymmetric). (a) p = q = 16, e _1 = 50, and t = 1.6 x 10 5 ; 
(b) (p,q) = (16,80), e _1 = 30, and t = 5.98 x 10 4 ; (c) 
(p,q) = (5,60), e- 1 = 25, and t = 4.89 x 10 5 . 
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This figure "figla.gif" is available in "gif" format from: 



http://arXiv.org/ps/patt-sol/9606002v2 
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http://arXiv.org/ps/patt-sol/9606002v2 
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